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This paper describes the development and analysis of finite-volume methods for the 
Landau-Lifshitz Navier-Stokes (LLNS) equations and related stochastic partial differential 
equations in fluid dynamics. The LLNS equations incorporate thermal fluctuations into 
macroscopic hydrodynamics by the addition of white-noise fluxes whose magnitudes are set 
by a fluctuation-dissipation relation. Originally derived for equilibrium fluctuations, the 
LLNS equations have also been shown to be accurate for non-equilibrium systems. Previous 
studies of numerical methods for the LLNS equations focused primarily on measuring vari- 
ances and correlations computed at equilibrium and for selected non-equilibrium flows. In 
this paper, we introduce a more systematic approach based on studying discrete equilibrium 
structure factors for a broad class of explicit linear finite- volume schemes. This new ap- 
proach provides a better characterization of the accuracy of a spatio-temporal discretization 
as a function of wavenumber and frequency, allowing us to distinguish between behavior at 
long wavelengths, where accuracy is a prime concern, and short wavelengths, where stability 
concerns are of greater importance. We use this analysis to develop a specialized third-order 
Runge Kutta scheme that minimizes the temporal integration error in the discrete structure 
factor at long wavelengths for the one-dimensional linearized LLNS equations. Together with 
a novel method for discretizing the stochastic stress tensor in dimension larger than one, our 
improved temporal integrator yields a scheme for the three-dimensional equations that sat- 
isfies a discrete fluctuation-dissipation balance for small time steps and is also sufficiently 
accurate even for time steps close to the stability limit. 
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I. INTRODUCTION 

Recently the fluid dynamics community has considered increasingly complex physical, chemical, 
and biological phenomena at the microscopic scale, including systems for which significant inter- 
actions occur across multiple scales. At a molecular scale, fluids are not deterministic; the state 
of the fluid is constantly changing and stochastic, even at thermodynamic equilibrium. As simula- 
tions of fluids push toward the microscale, these random thermal fluctuations play an increasingly 
important role in describing the state of the fluid, especially when investigating systems where the 
microscopic fluctuations drive a macroscopic phenomenon such as the evolution of instabilities, or 
where the thermal fluctuations drive the motion of suspended microscopic objects in complex flu- 
ids. Some examples in which spontaneous fluctuations can significantly affect the dynamics include 
the breakup of droplets in jets [U El |3], Brownian molecular motors [H El El [7] , Rayleigh-Bernard 
convection (both single species pj and mixtures [9], Kolmogorov fiows |10 | 111 ! [T2]. Ray leigh- Taylor 
mixing \13\ ITS] , combustion and explosive detonation |15| [T6] , and reaction fronts pT] . 

Numerical schemes based on a particle representation of a fiuid (e.g., molecular dynamics, di- 
rect simulation Monte Carlo [18]) inherently include spontaneous fiuctuations due to the irregular 
dynamics of the particles. However, by far the most common numerical schemes in computational 
fiuid dynamics are based on solving partial differential equations. To incorporate thermal fiuctu- 
ations into macroscopic hydrodynamics. Landau and Lifshitz introduced an extended form of the 
compressible Navier-Stokes equations obtained by adding white-noise stochastic fiux terms to the 
standard deterministic equations. While they were originally developed for equilibrium fiuctua- 
tions, specifically the Rayleigh and Brillouin spectral lines in light scattering, the validity of the 
Landau-Lifshitz Navier-Stokes (LLNS) equations for non-equilibrium systems has been assessed 
|19| and verified in molecular simulations |20| [21] [22] . The LLNS system is one of the more com- 
plex examples in a broad family of PDEs with stochastic fiuxes. Many members of this family 
arise from the LLNS equations in a variety of approximations (e.g., stochastic heat equation) while 
others are stochastic variants of well-known PDEs, such as the stochastic Burger's equation [23], 
which can be derived from the continuum limit of an asymmetric excluded random walk. 

Several numerical approaches for fiuctuating hydrodynamics have been proposed. The earli- 
est work by Garcia et al. [23] developed a simple scheme for the stochastic heat equation and 
the linearized one-dimensional LLNS equations. Ladd et al. have included stress fiuctuations in 
(isothermal) Lattice-Boltzmann methods for some time [25], and recently a better theoretical foun- 
dation has been established j26| I27j. Moseler and Landman fT] included the stochastic stress tensor 
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of the LLNS equations in the lubrication equations and obtain good agreement with their molecular 
dynamics simulation in modeling the breakup of nano-jets. Sharma and Patankar [28] developed 
a fluid-structure coupling between a fluctuating incompressible solver and suspended Brownian 
particles. Coveney, De Fabritiis, Delgado-Buscalioni and co-workers have also used the isothermal 
LLNS equations in a hybrid scheme, coupling a continuum fluctuating solver to a Molecular Dy- 
namics simulation of a liquid \29\ [501 EI] • Atzberger and collaborators [32 have developed a version 
of the immersed boundary method that includes fluctuations in a pseudo-spectral method for the 
incompressible Navier-Stokes equations. Voulgarakis and Chu [33] developed a staggered scheme 
for the isothermal LLNS equations as part of a multiscale method for biological applications, and 
a similar staggered scheme was also described in Ref. |34j . 

Recently, Bell et al. [35] introduced a centered scheme for the LLNS equations based on interpo- 
lation schemes designed to preserve fluctuations combined with a third-order Runge-Kutta (RK3) 
temporal integrator. In that work, the principal diagnostic used for evaluation of the numerical 
method was the accuracy of the local (cell) variance and spatial (cell-to-cell) correlation structure 
for equilibrium and selected non-equilibrium scenarios (e.g., constant temperature gradient). The 
metric established by those types of tests is, in some sense, simultaneously too crude and too 
demanding. It is too crude in the sense that it provides only limited information from detailed 
simulations that cannot be directly linked to specific properties of the scheme. On the other hand, 
such criteria are too demanding in the sense that they place requirements on the discretization in- 
tegrated over all wavelengths, requiring that the method perform well at high wavenumbers where 
a deterministic PDE solver performs poorly. Furthermore, although Bell et al. [35^ demonstrate 
that RK3 is an effective algorithm, compared with other explicit schemes for the compressible 
Navier-Stokes equations, the general development of schemes for the LLNS equations has been 
mostly trial-and-error. 

Here, our goal is to establish a more rational basis for the analysis and development of explicit 
finite volume scheme for SPDEs with a stochastic flux. The approach is based on analysis of 
the structure factor (equilibrium fluctuation spectrum) of the discrete system. The structure 
factor is, in essence, the stationary spatio-temporal correlations of hydrodynamic fluctuations as a 
function of spatial wavenumber and temporal frequency; the static structure factor is the integral 
over frequency (i.e., the spatial spectrum). By analyzing the structure factor for a numerical 
scheme, we are able to develop notions of accuracy for a given discretization at long wavelengths. 
Furthermore, in many cases the theoretical analysis for the structure factor is tractable (with the 
aid of symbolic manipulators) allowing us to determine optimal coefficients for a given numerical 
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scheme. We perform this optimization as a two-step procedure. First, a spatial discretization is 
developed that satisfies a discrete form of the fluctuation-dissipation balance condition. Then, a 
stable temporal integrator is proposed and the covariances of the random numbers are chosen so as 
to maximize the order of temporal accuracy of the small-wavenumber static structure factor. We 
focus primarily on explicit schemes for solving the LLNS equations because even at the scales where 
thermal fluctuations are important, the limitation on time step imposed by stability is primarily 
due to the hyperbolic terms. That is, when the cell size is comparable to the length scale for 
molecular transport (e.g., mean free path in a dilute gas) the time step for these compressible 
hydrodynamic equations is limited by the acoustic CFL condition. At even smaller length scales 
the viscous terms further limit the time step yet the validity of a continuum representation for the 
fluid starts to break down at those atomic scales. 

The paper is divided into roughly two parts: The first half (sections [ll IV) defines notation. 



develops the formalism, and derives the expressions for analyzing a general class of linear stochastic 
PDEs from the LLNS family of equations. The main result in the first half, how to evaluate the 



structure factor for a numerical scheme, appears in section IIIB The second half applies this 



analysis to systems of increasing complexity, starting with the stochastic heat equation (section 



VA), followed by the LLNS system in one dimension (section VI) and three dimensions (section 



VII). The paper closes with a summary and concluding remarks. 



II. LANDAU-LIFSHITZ NAVIER-STOKES EQUATIONS 

We consider the accuracy of explicit finite-volume methods for solving the Landau-Lifshitz 
Navier-Stokes (LLNS) system of stochastic partial differential equations (SPDEs) in d dimensions, 
given in conservative form by 

dtU = -V-[F{U)-Z{U,r,t)], (1) 

r 1 

where U {r,t) = e is a vector of conserved variables that are a function of the spatial 

position r and time t. The conserved variables are the densities of mass p, momentum j = pv, and 
energy e = e{p,T) + \pv^, expressed in terms of the primitive variables, mass density p, velocity 
V and temperature T; here e is the internal energy density. The deterministic flux is taken from 
the traditional compressible Navier-Stokes-Fourier equations and can be split into hyperbolic and 
diffusive fluxes: 

F{U)=Fh{U) + Fd{U), 
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where 



H 



pv 

pvv^ + PI 
{e + P)v 



and Ff) 




a 

(T v + $, 



P = P{p, T) is the pressure, the viscous stress tensor \s ct = 2ri \{^v + Vi» 



d 



for d > 2 



(we have assumed zero bulk viscosity) and a = rjVx for d = 1, and the heat flux is ^ = /uVT. We 
denote the adjoint (conjugate transpose) of a matrix or linear operator M with M* = M . As 
postulated by Landau-Lifshitz |19| [36] , the stochastic flux 



is composed of the stochastic stress tensor Xl and stochastic heat flux vector S, assumed to be 
mutually uncorrelated random Gaussian fields with a covariance 
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(I](r, t)I]*(r', t')) =Cj:5{t - t')S{r - r'), where c{g« = '^vksT (^S^kSji + 5^5, 
(H(r,t)H^(r',t')> =Cs5{t - t')5{r - r'), where cff = 2JlkBT^5^J 



d 



f 



(2) 



In the LLNS system, the hyperbolic or advective fluxes are responsible for transporting the 
conserved quantities at the speed of sound or fluid velocity, without dissipation. On the other hand, 
the diffusive or dissipative fluxes are the ones responsible for damping the thermal fluctuations 
generated by the stochastic or fluctuating fluxes. At equilibrium a steady state is reached in which 
a fluctuation-dissipation balance condition is satisfied. 

In the original formulation. Landau and Lifshitz only considered adding stochastic fluxes to the 
linearized Navier-Stokes equations, which leads to a well-defined system of SPDEs whose equilib- 
rium solutions are random Gaussian fields. Derivations of the equations of fluctuating hydrody- 
namics through careful asymptotic expansions of the underlying microscopic (particle) dynamics 
give equations for the Gaussian fluctuations around the solution to the usual deterministic Navier- 
Stokes equations [37j, in the spirit of the Central Limit Theorem. Therefore, numerical solutions 
should, in principle, consist of two steps: First solving the nonlinear deterministic equations for 
the mean solution, and then solving the linearized equations for the fluctuations around the mean. 
If the fluctuations are small perturbations, it makes sense numerically to try to combine these two 
steps into one and simply consider non-linear equations with added thermal fluctuations. There 
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is also hope that this might capture effects not captured in the two-system approach, such as 
the effect of fluctuations on the very long-time dynamics of the mean (e.g., shock drift [35]) or 
hydrodynamic instabilities [HISIE!- 

The linearized equations of fluctuating hydrodynamics can be given a well-defined interpretation 
with the use of generalized functions or distributions [38]. However, the non-linear fluctuating 
hydrodynamic equations ([T]) must be treated with some care since they have not been derived 
from first-principles [19 and are in fact mathematically ill-defined due to the high irregularity of 
white-noise fluctuating stresses [39] . More specifically, because the solution of these equations is 
itself a distribution the interpretation of the nonlinear terms requires giving a precise meaning 
to products of distributions, which cannot be defined in general and requires introducing some 
sort of regularization. Although written formally as an SPDE, the LLNS equations are usually 
interpreted in a finite volume context, where the issues of regularity, at first sight, disappear. 
However, in finite volume form the level of fluctuations becomes increasingly large as the volume 
shrinks and the non-linear terms diverge leading to an "ultraviolet catastrophe" of the kind familiar 
in other fields of physics. Furthermore, because the noise terms are Gaussian, it is possible for rare 
events to push the system to states that are not thermodynamically valid such as negative T or p. 
For that reason, we will focus on the linearized LLNS equations, which can be given a well-defined 
interpretation. Since the fluctuations are expected to be a small perturbation of the deterministic 
solution, the nonlinear equations should behave similarly to the linearized equations anyway, at 
least near equilibrium for sufficiently large cells. 

To simplify the exposition we assume the fluid to be a mono- atomic ideal gas; the generalization 
of the results for an arbitrary fluid is tedious but straightforward. For an ideal gas the equation 
of state may be written as -P = p{kBT/m) = pc^, where c is the isothermal speed of sound. 
The internal energy density is e = pc^T, where is the heat capacity at constant volume, which 
may be written as = dfkB/2m where dj is the number of degrees of freedom of the molecules 
(for monoatomic gases there are dj = d translational degrees of freedom), and Cp = {1 + 2/df)c^ 
is the heat capacity at constant pressure. For analytical calculations it is convenient to convert 
the LLNS system from conserved variables to primitive variables, since the primitive variables are 
uncorrelated at equilibrium and the equations ([T]) simplify considerably. 
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where DfO = dfO + • V (□) denotes the famihar advective derivative. Note that in the fully non- 
linear numerical implementation, however, we continue to use the conserved variables to ensure 
that the physical conservation laws are strictly obeyed. 

Linearizing ([3| around a reference uniform equilibrium state p = po + 5p, v = vq + 5v, T = 
To + (5T, and dropping the deltas for notational simplicity. 
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we obtain the linearized LLNS system for the equilibrium thermal fluctuations, 
dtU = -V • [FU -Z] = -V- [FhU + FdVU - Z] , 

where 



(4) 



FhU 



PqV + pvo 
[cIpq^P + cITq'T)I + vov^ 
clc-^v + Tvo 



and FdVU 





Po^iloVv 



^(V^; + Vv^) - ^^I. Here Z{r, t) is a 



and V denotes a symmetrized traceless gradient, Vt? 
random Gaussian field with a covariance 

{Z{r, t)Z\r\ t')) = CzS{t - t')6{r - r' 

where the covariance matrix is block diagonal, 


Cz= 

Po^c-^Ch 

and Cs and Cs are given in Eq. (|2|. Equation Q is a system of linear SPDEs with additive noise 
that can be analyzed within a general framework, as we develop next. We note that the stochastic 
"forcing" in Q is essentially a divergence of white noise, modeling conservative intrinsic (thermal) 
fluctuations ^37j , rather than the more common external fluctuations modeled through white noise 
forcing [30]. 

The next two sections develop the tools for analyzing explicit finite volume schemes for lin- 
earized SPDEs, such as the LLNS system, specifically how to predict the equilibrium spectrum of 
the fluctuations (i.e., structure factor) from the spatial and temporal discretization used by the 



numerical algorithm. These analysis tools are demonstrated for simple examples in Section VA 



and applied to the LLNS system in Sections VI and VII 
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III. EXPLICIT METHODS FOR LINEAR STOCHASTIC PARTIAL DIFFERENTIAL 

EQUATIONS 

In this section, we develop an approach for analyzing the behavior of explicit discretizations for 
a broad class of SPDEs, motivated by the linearized form of the LLNS equations. In particular, 
we consider a general linear SPDE for the stochastic field Vl{r, t) = tl{t) of the form 

dU{t) = CU{t)dt + KdB{t), (5) 

with periodic boundary conditions on the torus r G V = [0,ff]'^, where C (the generator) and /C 
(the filter) are time-independent linear operators, and B is a cylindrical Wiener process (Brownian 
sheet), and the initial condition at i = is IAq. As common in the physics literature, we will abuse 
notation and often write 

dtU = CU + /CW, 

where W = d3{t)/dt is spatio-temporal white noise, i.e., a random Gaussian field with zero mean 
and covariance 

(W(r, t)yV*(r', t')) = 5{t - t')5{r - r'). (6) 

The so-called mild solution [38j of ([5]) is a generalized process 

U{t) = e^^Uo + /* e^*-'^^KdB{s), (7) 
Jo 

where the integral denotes a stochastic convolution. If the operator C is dissipative, that is, 
e*^Wo ~^ for all IAq, then at long times t' the solution to (5) is a Gaussian process with mean 

t—*co 

zero and covariance 

Cu{t) = {U{t')U*{t' + t))= / e-'^JCK.*e^'-'^^*ds, t > 0. (8) 

This means that ([5]) has a unique invariant measure (equilibrium or stationary distribution) that 
is Gaussian with mean zero and covariance given in Eq. ([8|. 

In general, the field Vl{r, t) is only a generalized function of the spatial coordinate r and cannot 
be evaluated pointwise. For the cases we will consider here, specifically, translationally-invariant 
problems where C and /C are differential operators, this difficulty can be avoided by transforming 
^ to Fourier space via the Fourier series transform 

U{r,t)=^e''^-^a{k,t) (9) 



fcGV 



1 



V 



U{k,t) = ^l e-''^-''U{r,t)dr, (10) 



rev 
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where = |V| = H'^ is the volume of the system, and each wavevector k = k{K) is expressed in 
terms of the integer wave index n G Z'^, giving the set of discrete wavevectors 



In Fourier space, the SPDE ^ becomes an infinite system of uncoupled stochastic ordinary dif- 
ferential equations (SODEs), 

(£l{t) = CU{t)dt + kdB{t), (11) 



one SODE for each fc G V . The invariant distribution of (11) is a zero-mean Gaussian random 

process, characterized fully by the covariance obtained from the spatial Fourier transform of ([s]), 

/ \ 1 

S{k,t) = V {U{k,t')U{k,t' + t)) = — e''^^S{k,Lo)du;, (12) 
\ / 2tt 

where the dynamic structure factor (space-time spectrum) is 

S{k,iv) = V (u{k,uj)U*{k,uj)'^ = (2-iu?j ^ [kK^ [jC + iu?j \ (13) 

which follows directly from the space-time {k, ui) Fourier transform of the SPDE ^ . By integrating 
the dynamic spectrum over all frequencies to, one gets the static structure factor 

1 r°° 

S{k) =S{k,t = 0) = — S{k,Lu)dio, (14) 

which is the spatial spectrum of an equilibrium snapshot of the fluctuating field and is the Fourier 
equivalent of Cu{t = 0). Note that the static structure factor of spatial white noise (a snapshot of 
W) is unity independent of the wavevector, «Sw;(^) = V (VV(fc, t)W*{k, t)) = I. 



A. Discretization 



For the types of equations we will consider in this paper, the invariant measure is spatially white, 
specifically, S{k) is diagonal and independent of k. The associated fluctuating field lA cannot be 
evaluated pointwise, therefore, it is more natural to use finite-volume cell averages, denoted here by 
U. In the deterministic setting, for uniform periodic grids there is no important difference between 
finite-volume and finite-difference methods. Our general approach can likely be extended also to 
analysis of stochastic finite-element discretizations, however, such methods have yet to be developed 
for the LLNS equations and here we focus on finite-volume methods. For notational simplicity, we 
will discuss problems in one spatial dimension (d = 1), with (mostly) obvious generalizations to 
higher dimensions. 
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Uj{t) = — / U{x,t)dx. (15) 



Space is discretized into Nc identical cells of length Ax = H/Nc, and the value Uj stored in 
cell 1 < i < A'^c is the average of the corresponding variable over the cell 

'{j-l)Ax 

Time is discretized with a time step At, approximating cell averages of I4{x, t) pointwise in time 
with U^ = {U1,...,U%^}, 

U] « UjinAt), 

where re > enumerates the time steps. The white noise W(x, t) cannot be evaluated pointwise 
in either space or time and is discretized using a spatio-temporal average 



1 />(n+l)At pAx 

^» = 7wu/ / mx,t)dxdt, (16) 

which is a normal random variable with zero mean and variance (AxAt)^^, independent between 
different cells and time steps. Note that for certain types of equations the dynamic structure factor 
may be white in frequency as well. In this case, a pointwise-in-time discretization is not appropriate 



and one can instead use a spatio-temporal average as done for white noise in (16). 

We will study the accuracy of explicit linear finite- volume schemes for solving the SPDE 
Rather generally, such methods are specified by a linear recursion of the form 

= (I + LAt) [/" + (17) 

where L and K are consistent stencil discretizations of the continuum differential operators C and 
K, (note that L and K may involve powers of At in general). Here 

VT" = (AxAt)5 (18) 

is a vector of standard normal variables with mean zero and variance one. 

Without the random forcing, the deterministic equation tit = CXi and the associated discretiza- 
tion can be studied using classical tools and notions of stability, consistency, and convergence. 
Under the assumption that the discrete generator L is dissipative, the initial condition will be 
damped and the equilibrium solution will simply be a constant. The addition of the random forc- 
ing, however, leads to a non-trivial invariant measure (equilibrium distribution) of [/" determined 
by an interplay between the (discretized) fluctuations and dissipation. Because of the dissipative 
nature of the generator, any memory of the initial condition will eventually disappear and the 
long time dynamics is guaranteed to follow an ergodic trajectory that samples the unique invariant 
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measure. In order to characterize the accuracy of the stochastic integrator, we will analyze how 
well the discrete invariant measure (equilibrium distribution) reproduces the invariant measure of 
the continuum SPDE (this is a form of weak convergence). Note that due to ergodicity, ensemble 
averages can either be computed by averaging the power spectrum of the fields over multiple sam- 
ples or averaging over time (after sufficiently many initial equilibration steps). In the theory we 
will consider the limit n ^ oo and then average over different realizations of the noise W to obtain 
the discrete structure factors. In numerical calculations, we perform temporal averaging. 



Regardless of the details of the iteration (17), VF" will always be a Gaussian random vector 



generated anew at each step n using a random number generator. The discretized field U"^ is 
therefore a linear combination of Gaussian variates and it is therefore a Gaussian vector-valued 
stochastic process. In particular, the invariant measure (equilibrium distribution) of U"" is fully 
characterized by the covariance 

which we would like to compare to the covariance of the continuum Gaussian field Cu {t = nAi) 
given by (|8|. This comparison is best done in the Fourier domain by using the spatial discrete 
Fourier transform, defined for a spatially-discrete field U [for example, U = ?7" or U = U{t)] via 

U, = ^ Uke'^"^" (20) 
=V ^ f^.+ie-'^'^'Ax, (21) 

j=0 

where we have denoted the discrete dimensionless wavenumber AA; = kAx = 2tth/Nc, and the 
wave index is now limited to the first Nc values, 

Vd = {k = 2ttk/H I < K < Nc} C IC. 

Since the fields are real-valued, there is a redundancy in the Fourier coefficients Uk because of 
the Hermitian symmetry between k and Nc — n (essentially, the second half of the wave indices 
correspond to negative fc), and thus we will only consider < k < \_Nc/2\ , giving a (Nyquist) cutoff 
wavenumber kmax ~ vr/Ax. 

^ n 

What we would like to compare is the Fourier coefficients of the numerical approximation, Uj^ , 
with the Fourier coefficients of the continuum solution, tik{t = nAt). The invariant measure of 
has zero mean and is characterized by the covariance obtained from the spatial Fourier transform 
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of (191 



V lim (tJ 



k 




Ns+n 



) 



•k 



) 



(22) 



From the definition of the discrete Fourier transform it follows that for small Ak, i.e., smooth 
Fourier basis functions on the scale of the discrete grid, TJk{t) converges to the Fourier coefficient 
Vl{k,t = nAt) of the continuum field. Therefore, Sk,n is the discrete equivalent (numerical approx- 
imation) to the continuum structure factor S{k,t = nAt). We define a discrete approximation to 
be weakly consistent if 



for any chosen k £ V and t. This means that, given a sufficiently fine discretization, the numerical 
scheme can accurately reproduce the structure factor for a desired wave index and time lag. An 
alternative view is that a convergent scheme reproduces the slow (compared to At) and large- 
scale (compared to Ax) fiuctuations, that is, it accurately reproduces the dynamic structure factor 
S{k, uj) for small Ak = kAx and Alu = uoAt. Our goal here is to quantify this for several numerical 
methods for solving stochastic conservation laws and optimize the numerical schemes by tuning 
parameters to obtain the best possible approximation to <S(fc, lo) for small k and cj. 
Much of our analysis will be focused on the discrete static structure factor 



Note that for a spatially- white field Vl{x), the finite- volume averages Uj are independent Gaussian 
variates with mean zero and variance Ax^^, and the discrete Fourier coefficients Uk are independent 
Gaussian variates with mean zero and variance V~^. As a measure of the accuracy of numerical 
schemes for solving Eq. (|5]), we will compare the discrete static structure factors with the 
continuum prediction S{k), for all of the discrete wavenumbers (i.e., pointwise in Fourier space). 
It is expected that any numerical scheme will produce some artifacts at the largest wavenumbers 
because of the strong corrections due to the discretization; however, small wavenumbers ought to 
have much smaller errors because they evolve over time scales and length scales much larger than 
the discretization step sizes. Specifically, we propose to look at the series expansions 



and optimize the numerical schemes by maximizing the powers pi and p2- Next we describe 
the general formalism used to obtain explicit expressions for the discrete structure factors Sk for 
a general explicit method, and then illustrate the formalism on some simple examples, before 
attacking the more complex equations of fiuctuating hydrodynamics. 



lim S, 



'k,n=lt/At\ = «5 {k,t) 




Sk-S{k) = 0{AtP'kP^) 
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B. Analysis of Linear Explicit Methods 

Regardless of the details of a particular scheme and the particular linear SPDE being solved, at 
the end of the timestep a typical explicit scheme makes a linear combination of the values in the 
neighboring cells and random variates to produce an updated value, 

Aj=wr) Aj=tos 

U]+' = U^+ Yl *A,f/^A,+ E *A,^^A„ (23) 

Aj=—W£) Aj=—ws 

where wd and ws are the deterministic and stochastic stencil widths. The particular forms of the 
matrices of coefficients $ and ^ depend on the scheme, and will involve powers of At and Ax. 
Here we assume that for each n the random increment is an independent vector of Ns normal 
variates with covariance Cw = )*) constant for all of the cells j and thus wavenumbers, 

where Ng is the total number of random numbers utilized per cell per stage. Computer algebra 



systems can be used to obtain explicit formulas for the matrices in (23); we have made extensive 
use of Maple for the calculations presented in this paper. 



Assuming a translation invariant scheme, the iteration (23) can easily be converted from real 
space to an iteration in Fourier space, 

Aj=wr> Aj=ws 

UT^ = UI+ ^AjUlexp{iAjAk)+ Y ^^AjWlexpizAjAk), (24) 

Aj=-wo Aj=-ws 

where different wavenumbers are not coupled to each other. In general, any linear explicit method 
can be represented in Fourier space as a recursion of the form 

C/r' = MkUl + NkWl, (25) 



where the explicit form of the matrices Mk and Nk depend on the particular scheme and typically 

'w 



contain various powers of sinA/c, cosAA;, and At, and C^^ = (Wj^ ( VF;. ) ) = ^Cw By 



iterating this recurrence relation, we can easily obtain (assuming Uj. = 0) 

n 

V — ^ , ■ — n—l 
1=0 

from which we can calculate 

n— 1 n—l 

SI = V{ (Ul) (Uiy) = J2 (Mk)' {AxN.CwND {MlY = ^ (M,)^ C {MlY . 

1=0 1=0 

In order to calculate this sum explicitly, we will use the following identity 

MkSlMl -Sl = {MkT C {Mlf - C 
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to obtain a linear system for the entries of the matrix S^. If the deterministic method is stable, 
which means that all eigenvalues of the matrix Mj. are below unity for all wavenumbers, then in 
the limit n — > oo the first term on the right hand side will vanish, to give 

MkSkMl -Sk = -AxNkCwNl. (26) 

This is a linear system of equations for the equilibrium static structure factor produced by a given 
scheme, where the number of unknowns is equal to the square of the number of variables (field 
components). By simply deleting the subscripts k one obtains a more general but much larger 
linear system [5l] for the real space equilibrium covariance of a snapshot of the discrete field 



MCuM* -Cu = -AxNC'^^^N* 



where C^^^^ = (W)*) is the covariance matrix of the random increments. Note that this 
relation continues to hold even for schemes that are not translation invariant such as generalizations 
to non-periodic boundary conditions; however, the number of unknowns is now the square of the 
total number of degrees of freedom so that explicit solutions will in general not be possible. Based 
on standard wisdom for deterministic schemes, it is expected that schemes that perform well 
under periodic boundary conditions will also perform well in the presence of boundaries when the 
discretization is suitably modified only near the boundaries. 

A similar approach to the one illustrated above for the static structure factor can be used to 
evaluate the discrete dynamic structure factor 

from the time-discrete Fourier transform 

where Auj = uAt, and the frequency is less than the Nyquist cutoff, to < LOmax = vr/At. The 
calculation yields 

Sk^u^ = [I- exp {-iAuj) Mfc]-^ {AxAtNkCwND [I - exp {iAu) M^]"^ . (27) 



Equation ( |27| ) can be seen as discretized forms of the continuum version (13 1 in the limits Ak 0, 
At (the corresponding correlations in the time-domain are given in Ref. [41]). 



Equations (26 ) and (27) are the main result of this section and we have used it to obtain explicit 
expressions for and S^^uj for several equations and schemes. In the next sections we will illustrate 
the above formalism for several simple examples of stochastic conservation laws. 
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1. Discrete Fluctuation-Dissipation Balance 

Let us first consider the static structure factors for very small time steps. In the limit Ai — > 0, 
temporal terms of order two or more can be ignored so that all time-integration methods behave 



like an explicit first-order Euler iteration as in (17 1, 



uT'={l + ^^)OUsj^W,. (28) 

where = L (At = 0) can be thought of as the spatial discretization of the generator C, and 



1^(0) = K (At = 0) is the spatial discretization of the filtering operator /C. Comparing to (25) we 



can directly identify = I + AtL^^ ^ and = J ^k\, ^ and substitute these into Eq. (26). 



Keeping only terms of order At on both sides we obtain the condition 



t^>i^^sf{ii:'x=^'c^isf] . (29) 



where S^u^ = limAt^o (see also a related real-space derivation using Ito's calculus in Ref. 



as well as Section VIII in Ref. Pl]). It can be shown that if l\. ^ is definite, Eq. (29) has a 



unique solution. Assuming that W is as given in Eq. (18), i.e., that Cw = and that the spatial 



discretizations of the generator and filter operators satisfy a discrete fluctuation- dissipation balance 

£r+(£i°Y=-srf5frV. (30) 



we see that S^u^ = J is the solution to equation (29), that is, at equilibrium the discrete fields are 



spatially-white. The discrete fluctuation-dissipation balance condition can also be written in real 
space, 

^(o) + Mo)y^^(o)/^(o)y^ (31) 



The condition (31 ) is the discrete equivalent of the continuum fluctuation-dissipation balance con- 
dition ^ 

C + C* = -KK\ (32) 

which ensures that S{k) = /, i.e., that the invariant measure of the SPDE is spatially- white. 
We observe that adding a skew adjoint component to C does not alter the fluctuation-dissipation 
balance above, as is the case with non-dissipative (advective) terms. Numerous equations [S^ mod- 
eling conservative thermal systems satisfy condition ( |32[ ), including the linearized LLNS equations 
(with some additional prefactors). In essence, the fluctuations injected at all scales by the spatially 
white forcing W are filtered by /C and then dissipated by C at just equal rates. 
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IV. LINEAR STOCHASTIC CONSERVATION LAWS 

The remainder of this paper is devoted to the study of the accuracy of finite-volume methods 
for solving linear stochastic PDEs in conservation form, 

dtU = -V ■ [{AU - CVU) - EW] , (33) 

where A, C and E are constants, and W is Gaussian spatio-temporal white noise. The white 
noise forcing and its divergence here need to be interpreted in the (weak) sense of distributions 
since they lack the regularity required for the classical definitions. The linearization of the LLNS 
equations ([T| leads to a system of the form (33), as do a number of other classical PDEs [37 , such 



as the stochastic advection- diffusion equation 

dfT = -a-VT + iiV^T + y%V ■ W, (34) 

where T(r,t) = U{r,t) is a scalar stochastic field, A = a is the advective velocity, C = fil, fi > 
is the diffusion coefficient, and E = y/2fj,I. The simplest case is the stochastic heat equation, 
obtained by taking a = 0. 

A key feature the type of system considered here is that the noise is intrinsic to the system and 
appears in the flux as opposed to commonly treated systems that include an external stochastic 
forcing term, such as the form of a stochastic heat equation considered in Ref. [40 . Since white 
noise is more regular than the spatial derivative of white noise, external noise leads to more regular 
equilibrium fields (e.g., continuous functions in one dimension). Intrinsic noise, on the other hand. 



leads to very irregular equilibrium fields. Notationally, it is convenient to write (33) as 



dtU = -T> {AU - CgU - EW) , (35) 

defining the divergence "D = V- and gradient ^ = V operators, T)* = —Q. In the types of 
equations that appear in hydrodynamics, such as the LLNS equations, the operator T^A is skew- 
adjoint, iT>AY = — X>A (hyperbolic or advective flux), C ^ (dissipative or diffusive flux), and 
EE* = 2C, i.e., E* = {2Cfl'^. Therefore, the generator £ = -X>A T)CQ = (X>A)* - T)CT>* 



and filter /C = T^E satisfy the fluctuation-dissipation balance condition (32) and the equilibrium 



distribution is spatially- white. Note that even though advection makes some of the eigenvalues of 



C complex, the generator is dissipative and (33) has a unique invariant measure because the real 



part of all of the eigenvalues of C is negative except for the unique zero eigenvalue. 

It is important to point out that discretizations of the continuum operators do not necessarily 



satisfy the discrete fluctuation-dissipation condition (31). One way to ensure the condition is 
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satisfied is to discretize the diffusive components of the generator Ld = DCG and the filter 
K = DE using a discrete divergence D and discrete gradient G so that the discrete fluctuation- 
dissipation balance condition -Ld + L*^ = —KK* holds. If, however, the discretization of the 
advective component of the generator La = —DA is not skew-adjoint, this can perturb the balance 



(30). Notably, various upwinding methods lead to discretizations that are not skew-adjoint. The 
correction to the structure factor S^*^^ = I + ^S^k'' ^° ^ non-zero IS.La = {La + -L^) /2 can 



easily be obtained from Eq. (29), and in one dimension the result is simply 

We will use centered differences for the advective generator in this work, which ensures a skew- 
adjoint La, and our focus will therefore be on satisfying the discrete fluctuation-dissipation balance 
for the diffusive terms. 



Finite- Volume Numerical Schemes 



We consider here rather general finite- volume methods for solving the linear SPDE (33) in one 
dimension, 



dM = -^[Hu)-z] = -^ 

ox ox 



A-cS- \ U-EW 



(37) 



dx ^ 

with periodic boundaries, where we have denoted the modified (potentially correlated) white noise 
flux with 2 = EW. As for classical finite-volume methods for the deterministic case, we start 
from the PDE and integrate the left and right hand sides over a given cell j over a given time step 
At, and use integration by parts to obtain the formally exact 

where the deterministic discrete fluxes F and dimensionless discrete stochastic fluxes Z are calcu- 
lated on the boundaries of the cells (points in one dimension, edges in two dimensions, and faces in 
three dimensions), indexed here with half-integers. These fluxes represent the total rate of trans- 



port through the interface between two cells over a given finite time interval At, and (38) is nothing 
more than a restatement of conservation. The classical interpretation of pointwise evaluation of 
the fluxes is not appropriate because white noise forcing lacks the regularity of classical smooth 
forcing and cannot be represented in a finite basis. Instead, just as we projected the fluctuating 
fields using finite-volume averaging, we ought to project the fluxes to a finite representation as well 
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through spatio-temporal averaging, as done in Eq. (16). For the purposes of our analysis, one can 
simply think of the discrete fluxes as an approximation that has the same spectral properties as 
the corresponding continuum Gaussian fields over the wavevectors and frequencies represented by 
the finite discretization. 

The goal of numerical methods is to approximate the fluxes as best as possible. In general, 
within each time step of a scheme there may be Ngt stages or substeps; for example, in the classic 
MacCormack method there is a predictor and a corrector stage {N^t = 2), and in the three- 
stage Runge-Kutta method of Williams et al. [35] there are three stages {Ngt = 3). Each stage 



< s < Nst is of the conservative form ( 38 ) , 



-1 



^ - ^ " Ax [^^H " ^.-U + Ax37^ V " 

s'=0 



where the a's are some coefficients, ^i^^ — 1) ^^'^ each of the stage fluxes are partial approx- 

imations of the continuum flux. For the stochastic integrators we discuss here, the deterministic 
fluxes are calculated the same way as they would be in the corresponding deterministic scheme. In 
general, the stochastic fluxes Z-,i can be expressed in terms of independent unit normal variates 
• , 1 that are sampled using a random number generator. The stochastic fluxes in each stage may 
be the same, may be completely independent, or they may have non-trivial correlations between 
stages. 



Note that it is possible to avoid non-integer indices by re-indexing the fluxes in Eq. (38) and 



writing it in a form consistent with ( 23 ) , 

However, when considering the order of accuracy of the stencils and also fluctuation-dissipation 
balance in higher dimensions, it will become important to keep in mind that the fluxes are evaluated 
on the faces (edges or half-grid points) of the grid, and therefore we will keep the half-integer indices. 
Note that for face-centered values, such as fluxes, it is best to add a phase factor exp {iAk/2) in the 
definition of the Fourier transform, even though such pure phase shifts will not affect the correlation 
functions and structure factors. 

Before we analyze schemes for the complex LLNS equations, we present an illustrative expHcit 
calculation for the one-dimensional stochastic heat equation. 
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V. EXAMPLE: STOCHASTIC HEAT EQUATION 



We now illustrate the general formalism presented in Section IV for the simple case of an Euler 
and predictor-corrector scheme for solving the stochastic heat equation in one dimension, 



Vt = IJ-Vxx + y^Wx, 



(41) 



where v {x, t) =IA (x, t) is a scalar field and ^ is the mass or heat diffusion coefficient. The solution 
in the Fourier domain is trivial, giving 

2/iA;2 



u;2 + ^2^4 



, and S{k) = 1. 



(42) 



A. Static Structure Factor 



We first study a simple second-order spatial discretization of the dissipative fluxes 
combined with an Euler integration in time, to give a simple numerical method for solving the 



SPDE (41) 



(43) 



where u = U and the W^s are independent unit normal random numbers with zero mean generated 



anew at every time step (here Ns = Ngt = !)• From (43), we can extract the recursion coefficients 



appearing in (25 ), 



Mfc = 1 + Pie 



iAk 



2 + = l + 2/3(cosAA;-l), 



where 



AtV2 



denotes a dimensionless diffusive time step (ratio of the time step to the diffusive CFL limit). 



Together with Cw = l, (26) becomes a scalar equation for the discrete structure factor. 



{MkM^ - l)Sk = -AxNkN, 



ki 
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with dimensionless solution 



4/?(l — cosAfc) ^, ^, -,M-i 

Sk= (^_^2) = [l + /3(cosAfc-l)] 



(44) 



The time-dependent result can also easily be derived from (27) 



= f 1 - e-*/^) Sk, where i = nAt 



where t ^ = 4/i(cosAA; — 1) /Ax^ ~ 2/iA;^ is the familiar relaxation time for wavenumber k, 
showing that the smallest wavenumbers take a long time to reach the equilibrium distribution. 
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Figure 1: An illustration of the discrete structure factor for the Euler (43 1 [c.f Eq. (44)] and predictor 



corrector (45) [c.f. Eq. (44l] schemes for the stochastic heat equation (41) 



Equation (44) is a vivid illustration of the typical result for schemes for stochastic transport 
equations based on finite difference stencils, also shown in Fig. [T] Firstly, we see that for small 
k we have that S'fc ~ 1 + fi^k? jl^ showing that the smallest wavenumbers are correctly handled 
by the discretization for any time step. Also, this shows that the error in the structure factor is 
of order /?, i.e., of order At, as expected for the Euler scheme, whose weak order of convergence 
is one for SODEs. Finally, it shows that the error grows quadratically with k (from symmetry 
arguments, only even powers will appear). By looking at the largest wavenumber, Akmax = t^, we 
see that Sk^^^ = (1 — 2/3)"^, from which we instantly see the CFL stability condition /? < f/2 , 
which guarantees that the structure factor is finite and positive for all < A; < vr. Furthermore, 
we see that for /? <C 1, the structure factor is approximately unity for all wavenumbers. That is, a 
sufficiently small step will indeed reproduce the proper equilibrium distribution. 
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By contrast, a two-stage predictor-corrector scheme for the diffusion equation, 

(45) 

achieves much higher accuracy, namely, a structure factor that deviates from unity by a higher 
order in both At and k, 

PC-IRNG: Sk-1- 

as illustrated in Fig. [T} We can also use different stochastic fluxes in the predictor and the corrector 
schemes (i.e., use Ng = 2 random numbers per cell per stage), with an added pre- factor of \/2 to 
compensate for the variance reduction of the averaging between the two stages. 



-1 =-1 + (-"-1 - + + 2 {w^;4' - W^r) (predictor) 



< + ti" + ^ - 2n" + n^+i + 2V^— — W^r' - ' 



,(n,C) Tr.(n,C) 



(corrector). 

(46) 



For the scheme (46) the analysis reveals an even greater spatio-temporal accuracy of the static 



structure factors, namely, third order temporal accuracy 

PC-2RNG: Sk^l + (5^Ak^/8. 

This illustrates the importance of the handling of the stochastic fluxes in multi-stage algorithms, 
as we will come back to shortly. Note that the analysis we presented here for explicit methods can 
easily be extended to implicit and semi-implicit schemes as well, as illustrated in Appendix [T] for 
the Crank-Nicolson method for the stochastic heat equation. 

Previous studies |29| |35] have measured the accuracy of numerical schemes through the variance 
of the fields in real space, which, by Parseval's theorem, is related to the integral of the structure 



factor over all wavenumbers. For the Euler scheme (43) for the stochastic heat equation this can 
be calculated analytically, 

al = {u]) - {ujf = Ax-i (1 - 2/3)-i/2 ^ ^ ^ 



showing first-order temporal accuracy (in the weak sense). For the predictor-corrector scheme (45 ), 



on the other hand, (o"^*^)^ « Ax ^ (l — 3/3^/2). It is important to note, however, that using the 
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variance as a measure of accuracy of stochastic real-space integrators is both too rough and also 
too stringent of a test. It does not give insights into how well the equipartition is satisfied for the 
different modes, and, at the same time, it requires that the structure factor be good even for the 
highest wavenumbers, which is unreasonable to ask from a finite-stencil scheme. 

For pseudo-spectral methods, as studied for the incompressible fluctuating Navier-Stokes equa- 
tion in Ref. |44| , one can modify the spectrum of the stochastic forcing so as to balance the 
numerical stencil artifacts, and one can also use an (exact) exponential temporal integrator in 
Fourier space to avoid the artifacts of time stepping. However, for finite-volume schemes, a more 
reasonable approach is to keep the stochastic fluxes uncorrelated between disjoint cells (which is ac- 
tually physical), and instead of looking at the variance, focus on the accuracy of the static structure 
factor for small wavenumbers. Specifically, basic schemes will typically have — 1 = O (AtA:^), 
while multi-step schemes will typically achieve Sfc — 1 = O (At^A;^) or higher temporal order, or 
even 5^ - 1 = O (At'^k^). 



B. Dynamic Structure Factor 

It is also constructive to study the full dynamic structure factor for a given numerical scheme, 
especially for small wavenumbers and low frequencies. This is significantly more involved in terms 
of analytical calculations and the results are analytically more complicated, especially for multi- 



stage methods and more complex equations. For the Euler scheme (43) the solution to Eq. (27) 
is 

g ^ 2X1X2 V^^ 

2Ai-2(l -cosAcu) + x?X2 V^A;4' 

where xi = 2(1 — cosAk)/Ak'^ and X2 = 1 + 2/3 (cos AA; — 1). This shows that the dynamic 
structure factor does not converge to the correct answer for all wavenumbers even in the limit 
At 0, namely 

hm Sk,u> = — 2"T7i:- (47) 



For small AA;, xi ^ 1 — Ak^/6, and the numerical result closely matches the theoretical result (42 1. 
However, for finite wavenumbers the effective diffusion coefficient is multiplied by a prefactor xi, 
which represents the spatial truncation error in the second-order approximation to the Laplacian. 
For all of the time-integration schemes for the stochastic heat equation discussed above, one can 



23 



reduce the discrete dynamic structure factor to a form 

c, 2xstochlJ'k'^ 

~ 2At-2 (1 - cos Aw) + xL^^^^ ' 

where Xstoch and Xdet depend on (3 and Ak and can be used to judge the accuracy of the scheme. 

In this paper we focus on the static structure factors in order to optimize the numerical schemes 

and then simply check numerically that they also produce reasonably-accurate results for the 

dynamic structure factors for small and intermediate wavenumbers and frequencies. 



C. Higher-Order Differencing 

Another interesting question is whether using a higher-order differencing formula for the viscous 



fluxes improves upon the second-order formula in the basic Euler scheme (43). For example, a 



standard fourth order in space finite difference yields the modified Euler scheme 



(48) 



Repeating the previous calculation shows that 



lim Sk = 6[7 -cosAkr\ (49) 

demonstrating that the fluctuation-dissipation theorem is not satisfied for this scheme at the dis- 
crete level even for infinitesimal time steps. This is because the spatial discretization operators in 



(48) do not satisfy the discrete fluctuation dissipation balance. 



In order to obtain higher-order divergence and Laplacian stencils that satisfy (30) we can start 
from a higher order divergence discretization D and then simply calculate the resulting discrete 
Laplacian L = —DD*. Here D should be a fourth-order (or higher) difference formula that 
combines four face-centered values, two on each side of a given cell, into an approximation to the 
derivative at the cell center. Conversely, D* combines the values from four cells, two on each side 
of a given face, into an approximation to the derivative at the face center. A standard fourth-order 
finite-difference stencil for D produces the higher-order Euler scheme, 

, At^/^ / 1 q 9 1 \ 



for which S^k, 1 + /3A/i; /2, which is the same leading-order error as the basic Euler scheme (43). 



On the other hand, the dynamic structure factor for small time steps is as in Eq. (47) but now 
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Xi = (1 — cos AA;)(13 — cos AA;)/ (72AA;^) ~ 1 — 3Afc^/320, which shows the higher spatial order of 
the scheme. 



Note that in (50) both the discretization of the Laplacian and of the gradient are of higher 



spatial order than in (43), however, the Laplacian operator is not of the highest order possible 
for the given stencil width. We will not use higher-order differencing for the diffusive fluxes in 
this work in order to avoid large Laplacian stencils like the one above. Rather, we will use the 
traditional second-order discretization and focus on the time integration of the resulting system. 

D. Handling of Advection 

The analysis we illustrated here for the stochastic heat equation can be directly applied to the 



scalar advection-diffusion equation (34) in one dimension, 

vt = -avx + fJ.Vxx + \l^y^x- (51) 

For example, a second-order centered difference discretization of the advective term —avx leads to 
the following explicit Euler scheme 

uf' = u^-1 {u]_,, - + /? (n^., - 2u^ + + V^|^ {wj^.^ -W^.^, (52) 

where the dimensionless advective CFL number is 

aAt 

a = —— = Pr, 
Ax 

and r = aAx/fi is the so-called cell Reynolds number and measures the relative importance of 
advective and diffusive terms at the grid scale. Note that this scheme is unconditionally unstable 
when /U = 0, specifically, the stability condition is a^/2 < (3 < 1/2. 



For the Euler method (52) the analysis yields a structure factor 



(l-rV4) 
l-ar/2 ' 2(1 -ar/2)^ 



showing that even the smallest wavenumbers have the wrong spectrum for a finite time step when 
|r| > 0, which is unacceptable in practice since it means that even the slowly-evolving large-scale 
fluctuations are not handled correctly. Adding an artificial diffusion A;U = ;u|r|/2tO;U leads to an 
improved leading order error. 
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It is well-known that adding such an artificial diffusion is equivalent to upwinding the advective 
term and leads to much improved stability for large r as well^. 

The second-order predictor-corrector time stepping scheme can be applied when advection is 
included as well. If |r| > the leading order errors are 

2 

PC-IRNG: - I - (l - y) 

3 

PC-2RNG: Sk - I - ^AA;^, (53) 

8 

showing that PC-2RNG gives a more accurate discrete structure factor than PC-IRNG even if 
advection is included as well. Note that the predictor-corrector method is unconditionally unstable 



when // = 0. In Section VI A we analyze a three-stage Runge-Kutta scheme that has a small leading 



order error in but is also stable when a < 1 even if /i = 0. 



VI. LLNS EQUATIONS IN ONE DIMENSION 



In this section, we will consider the linearized LLNS system Q for a mono-atomic ideal gas 
in one spatial dimension, that is, where symmetry dictates variability along only the x axis. As 
explained in the Introduction, focusing on an ideal gas simply fixes the values of certain coefficients 
and thus simplifies the algebra, without limiting the generality of our analysis. We will arbitrarily 
choose the number of degrees of freedom per particle to be = 1, even though in most cases of 
physical interest d/ = 3 is appropriate; this merely changes some of the constant coefficients and 
does not affect our discussion. Explicitly, the one-dimensional linearized LLNS equations are 



(54) 



where the covariance matrices of the stochastic fluxes are = 2riokBTQ and Cs = 2//ofcB7'o • In 



dtp 




pov + pvo 














a 


4po^P + Cq^o"^^ + ^ov 




Po^VoVx 


d 




dtv 


dx 


dx 


dx 


dtT _ 








Po (^v PoTx _ 







Fourier space the flux becomes 



^^0 Po 
cgc-i 



Tq 'c ; 



1„2 




[vo - ikpQ Vo) 



^ Note that for this particular type of upwinding the denominator in Eq. ( 36 1 vanishes identically and it can be 
shown that the correct solut 
discretizations of advection. 



shown that the correct solution is AS'^'^^ — 0, however, this is not necessarily true for other, higher order, upwind 
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which through Eqs. (13) and (14) (or, equivalently, Eq. (29)) gives static structure factors that 
are independent of k, 



S{k) 



PoCQ'^ksTo 
pQ^ksTo 









Po^c-^kBT^ 



(55) 



Therefore, the invariant distribution for the spatial fluctuating fields is white noise, uncorrelated 



among the different primitive variables, and with variances given in Eq. (55). This is in agreement 
with predictions of statistical mechanics, and how Landau and Lifshitz obtained the form of the 
stochastic fluxes. Note that in the incompressible limit, cq — > oo, the density fluctuations diminish, 
but the velocity and temperature fluctuations are independent of cq. 

In this section we will calculate the discrete structure factor for several finite-volume approxi- 



mations to ( 54 ) . From the diagonal elements of Sk we can directly obtain the non-dimensionalized 
static structure factors for the three primitive variables, for example, 

Sl'^ = 1 ^ iPkPl) , 
PoCq kslo 

which for a perfect scheme would be unity for all wavevectors. Similarly, the off-diagonal or cross 
elements, such as for example 

V 



{pkvD , 



{poCo'^kBTo) {pQ^kBTo) 

would all vanish for all wavevectors for a perfect scheme. Our goal will be to quantify the deviations 
from "perfect" for several methods, as a function of the discretization parameters Ax and At. 



A. Third-order Runge-Kutta (RK3) Scheme 

When designing numerical schemes to integrate the full LLNS system, it seems most appropri- 
ate to base the scheme on well-known robust deterministic methods, and modify the deterministic 
methods by simply adding a stochastic component to the fluxes, in addition to the usual determin- 
istic component. With such an approach, at least we can be confident that in the case of weak noise 
the solver will be robust and thus we will not compromise the fluid solver just to accommodate the 
fluctuations. 

A well-known approach to solving PDEs in conservation form 
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is to use the method of lines to decouple the spatial and temporal discretizations. We will focus 
on one dimension first for notational simplicity. In the method of lines, a finite-volume spatial 
discretization is applied to the obtain a system of stochastic differential equations for the discretized 
fields: 

dUi 



-Ax 



-1 



dt 



FniUj^i) - Fh{U^_i)] - Ax-' \Fd{V.^iU) - Fd{V^_iU)\ , (56) 



where U-,i are face-centered values of the fields that are calculated from the cell-centered values 
Uj, and V • , 1 is a cell-to- face discretization of the gradient operator. Any classical temporal 
integrator can be applied to the resulting system of SODEs. It is well known that the Euler and 
Heun (two-step second-order Runge-Kutta) methods are unconditionally unstable for hyperbolic 
equations. In Ref. [35 , an algorithm for the solution of the LLNS system of equations ([T| was 
proposed, which is based on the three-stage, low-storage TVD Runge-Kutta (RK3) scheme of 
Gottlieb and Shu |l6] . The RK3 scheme is the simplest TVD RK discretization for the deterministic 
compressible Navier-Stokes equations that is stable even in the inviscid limit, with the omission of 
slope-limiting. Here we adopt the same basic scheme and investigate optimal ways of evaluating 
the stochastic flux. 

In the RK3 scheme, the hyperbolic component of the face flux Fh is calculated by a cubic 
interpolation of U from the cell centers to the faces using an interpolation formula borrowed from 
the PPM method [47j, 

t/,.+ i = ^ {U, + J7,+i) - ^ (C/,_i + , (57) 

and then directly evaluating the hyperbolic flux from the interpolated values. In Refs. |35| SB] a 
modified interpolation is proposed that preserves variances; however, our analytical calculations 
indicate that this type of interpolation artificially increases the structure factor for intermediate 
wavenumbers in order to compensate for the errors at larger wavenumbers. Note that for the full 
non-linear equations, the conserved quantities are interpolated and then primitive face variables 
are calculated from those. For the linearized equations it does not matter and it is simpler to 
work exclusively with primitive variables. In the RK3 method, the diffusive components of the 
fluxes F £) are calculated using classical face-centered second-order centered stencils to evaluate the 
gradients of the fields at the cell faces. Stochastic fluxes Z-,i are also generated at the faces of 

-'"'"2 

the grid using a standard random number generator (RNG). These stochastic fluxes are generated 
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independently for velocity and temperature, and are zero for density, 





f.{RNG) 



2 

J^2 J 



where W 



(1/2) 



denotes a normal variate with zero mean and unit variance. 



For each stage of the RK3 scheme, a total cell increment is calculated as 



/^Uj{U,W) 



At 

Ax 



F^,^(U)-F^^{U) 



Each time step of the RK3 algorithm is composed of three stages 



=[/" + AC/j([/", Wi) (estimate at t = (n + l)At 



U 



n+l 



4^ 4 
3^ ^ 3 



(estimate at t = (n + 2)^^ ) 



J7'""5+Ai7,(J7"+3,iy3) 



(58) 



where for now we have not assumed anything about how the stochastic fluxes between different 
stages, Wi, W2 and W3, are related to each other. The relevant dimensionless parameters that 
measure the ratio of the time step to the CFL stability limits are 

coAi 



a 



Ax 
. r?oAt 

poAx"^ 

poAt 



a 
r 



1 a 



a 
P 



poCyAx^ Pr r 

where r = copoAx/rjQ is the cell Reynolds number and measures the relative importance of acoustic 
and viscous terms at the grid scale (we have assumed a low Mach number flow, i.e., \vq\ <^ cq), 
and Pr = r]oCv/po is the Prandtl number of the fluid. For low-density gases, r and p = rPr can be 
close to or smaller than one, however, for dense fluids sound dominates and r > 1 and p > 1 for all 
reasonable Ax (essentially. Ax > A, where A is the mean free path). In practice, in order to fully 
resolve viscous scales, one should keep both r and p reasonably small. 



B. Evaluation of the Stochastic Fluxes 



, a different stochastic flux is generated in each stage, that 
is, Wg = V^W^^^pj^, s = 1 ... 3. The additional prefactor \/2 is added because the averaging 



In the original RK3 algorithm 

As) „ 

RNG^ * 
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between the three stages reduces the variance of the overall stochastic flux. One can also use 

(s) 

different weights for each of the three stochastic fluxes, i.e., Ws = WgW j^j^q. Another option is 
to simply use the same stochastic flux W^j^l^Q in all three stages, that is, Ws = ^^^fG- A further 
option is to use the same random flux W^^j^Q in all three stages, but put in different weights in 
each stage, i.e., Wg = WsW^j^j^Q. Our goal is to find out which approach is optimal. For this 
purpose, we can generally assume that the three random fluxes are different, to obtain a total of 



six random numbers per cell per step, and use the formalism developed in Section III with A'^s = 6 
to express the structure factor in terms of the 6x6 covariance matrix of the random variates. 
This calculation is too tedious even for a computer algebra system, and we therefore first study 



the simple advection-diffusion equation (34) in order to gain some insight. 



1. Advection- Diffusion Equation 
The RK3 method can be directly applied to the scalar advection-diffusion equation in one 



dimension (51). Experience with deterministic solvers suggests that a numerical scheme that 
performs well on this type of model equation is likely to perform well on the full system ([T]) 
when viscous effects are fully resolved. Here we use the PPM-interpolation based discretization of 



the hyperbolic flux given in Eq. (57), which leads to a standard fourth-order centered difference 
approximation to the first derivative Vx [49^ (in Fourier space the relative error in the hyperbolic 
flux is of order 0{Ak^)), and thus justifies our choice for the interpolation. We discretize the 
gradient used in calculating the diffusive fluxes using the second-order centered difference 

V ■ , in — , 

-'^2 Ax 

which leads to the standard second-order centered difference approximation to the second derivative 
Vxx (the challenges with using the standard fourth-order centered difference approximation to Vxx 



are discussed in Section VC). The stencil widths in Eq. (23) are tu^i = 6 (three stages with 



stencil width two each) and ws = 4, and there are A^s = 3 random numbers per cell per step (one 



per stage), with a general 3x3 covariance matrix Cw. Equation (26) can then be solved to obtain 



the static structure factor for any wavenumber, however, these expressions are too complex to be 



useful for analysis. Instead, we perform an expansion of both sides of (26) for small k and thus 
focus on the behavior of the static structure factors for small wavenumbers and small time steps. 

As a first condition on Cu^, we have the weak consistency requirement Sfc=o = 1. With this 
condition satisfied, the method satisfies the discrete fluctuation-dissipation balance in the limit 
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At ^ since the discretization of the divergence is the negative adjoint of the discretization of 
the gradient. A second condition is obtained by equating the coefficient in front of the leading- 
order error term in S^, of order aAk"^, to zero; where the advective dimensionless CFL number is 
a = aAt/Ax. It turns out that this also makes the term of order aAk^ vanish. A third condition 
is obtained by equating the coefficient in front of the next-order error term of order a^Ak"^ to zero. 
Finally, a fourth condition equates the coefficient in front of a^AA;^ to zero. For this three-stage 
method, it is not possible to make the terms with higher powers of a vanish identically for any 
choice of Cw No additional conditions are obtained by looking at terms with powers of the 
diffusive CFL number (3 = /xAt/Ax^ since, as it turns out, the accuracy is always limited by the 
hyperbolic fluxes. 

The various ways of generating the stochastic fluxes can now be compared by investigating how 
many of these conditions are satisfied. It turns out that only the first condition is satisfied if we use 
a different independently-generated stochastic flux in each stage (one can satisfy one more condition 
by using different weights for the three independent stochastic fluxes). The second condition is 
satisfied if we use the same stochastic flux in all stages with a unit weight, i.e., Wg = WgW^^j^Q 
with wi = W2 = = 1. Armed with the freedom to put a different weight for this flux in each of 
the stages, we can satisfy the third condition as well if we use 



wi = W2 = = — , (59) 



which gives a structure factor 



Sk = l- ^,o?Ak^ - -^a^Ak'^ + h.o.t. 
24 6r^ 

If we are willing to increase the cost of each step and generate two random numbers per cell per 
step, we can satisfy the fourth condition as well. For this purpose, we look for a covariance matrix 
Cw that satisfies the four conditions and is also positive semi-definite and has a rank of two, i.e., 
has a smallest eigenvalue of zero. A solution to these equations gives the following method for 

evaluating the stochastic fluxes in the three stages 

^^3 =W^rL^ (60) 

where VF^g VK^^ are two independent random vectors that need to be generated and 
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stored during each RK3 step. This approach produces a structure factor 

T 24- -I- 

Sk = l- —a^Ak^ —a^Ak^ + h.o.t. 

24 288r 



We will refer to the RK3 scheme that uses one random flux per step and the weights in (59) as the 
RK3-1RNG scheme^ and to the RK3 scheme with two random fluxes per step as given in (60) as 
the RK3-2RNG scheme. 

It is important to point out that for the MacCormack method, which is equivalent to the 
Lax-Wendroff method for the advection-diffusion equation, the leading-order errors are of order 



aAk"^. This is much worse than for the stochastic heat equation (see Section V A) even though the 
MacCormack scheme is a predictor-corrector method. This is because of the low-order handling 
of advective fluxes used in the MacCormack method to stabilize the two-stage Runge-Kutta time 
integrator. 

C. Results for LLNS equations in One Dimension 

We can now theoretically study the behavior of the RK3-1RNG and RK3-2RNG schemes on the 



full linearized system (54), specializing to the case of zero background flow, = 0. As expected, 
we find that the behavior is very similar to the one observed for the advection-diffusion equation, 
in particular, the leading order terms have the same basic form. Specifically, the expansions of the 
diagonal and off-diagonal components of the structure factor Sk for the RK3-1RNG method are 

s\:^^sP^l + ^-^^l + e{a)Ak^ 



^i"'^^ ^i'-^a^Ak\ (61) 



where 



e(a) 



4 (3p + 2r) ■ 

These structure factors are shown in Fig. |2]for sample discretization parameters, along with the 
corresponding results for RK3-2RNG. We see from these expressions that as the speed of sound 
dominates the stability restrictions on the timestep more and more, namely, as p or r become larger 
and larger, a smaller a is required to reach the same level of accuracy, that is, a smaller timestep 
relative to the acoustic CFL stability limit is required. 
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Figure 2: An illustration of the discrete structure factor Sk for the LLNS equation for the RK3-1RNG (lines) 
and RK3-2RNG (same style of lines with added symbols) schemes, as calculated by numerical solution of 



(26 1 for an ideal one dimensional gas, for a = 0.5, /3 = 0.2 and (3t = 0.1. (Left) Diagonal (self) structure 
factors, which should ideally be identically unity. Also shown is the leading order error term 1 + e{a)Ak'^ 
(dotted line), which is the same for both schemes. (Right) Off-diagonal (cross) structure factors, which 
should ideally be identically zero. 



Similar results to Eqs. (61 ) hold also for the isothermal LLNS equations (in which the there is 



no energy equation), for which the calculations are simpler. For linearization around a constant 
background flow of speed vq = coMa, where Ma is the reference Mach number, the analysis for the 
isothermal LLNS equations shows that the error grows with the Mach number as 

Slf'^ ^ 1 + e{a) [1 + 6Ma2 + Ma^] A/fc^. 
VII. HIGHER DIMENSIONS 



Much of what we already described for one dimension applies directly to higher dimensions 
\35\ I48j . However, there is a peculiarity with the LLNS equations in three dimensions that does 
not appear in one dimension, and also does not appear for the scalar diffusion equation [42 . In 
one dimension the velocity component of the LLNS system of equations is essentially an advection- 
diffusion equation. In higher dimensions, however, there is an important difference, namely, the 
dissipation operator is a modified Laplacian Cm- By neglecting the hyperbolic coupling between 
velocity and the other variables in the linearized LLNS equations, we obtain the stochastic dijfusion 
equation 



^t = r]V- [C( Vi9)] + y2^V • \C^I'^W\ = r] {VCQ) i9 + y^T>C^/^W = rjCm^ + V^Wm, (62) 
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where C is the hnear operator that transforms the velocity gradient into a traceless symmetric 
stress tensor, 



C(Vi9) = 2 



2' '3 

and we have denoted the continuum velocity field with iD = 14 in order to distinguish from the 
discretized velocities v = U. Here we will focus on two-dimensional flows, i? = however, 
identical considerations apply to the fully three-dimensional case. 

If we arrange the components of the velocity gradient as a vector with four components 

1 T 



Vi9 



^x^xj ^x'^y^ ^y'^x^ ^y^y 



the linear operator C in ( 62 ) becomes the matrix 



C 



I -i 
110 
110 




(63) 



3 " " 3 

which is not diagonal. This means that the components of the stochastic stress C^^'^W would need 
to have non-trivial correlations between the x fluxes for Vx and y fluxes for Vy, as well as between 
the X fluxes for Vy and y fluxes for Vx- These correlations essentially amount to the requirement 
that the stochastic stress be a traceless symmetric tensor, at least at the level of its covariance 
matrix. Numerically, one generates independent random variates for the upper triangular portion 
of the stochastic stress tensor for each cell, then makes the tensor traceless and symmetric |19] . 
Note that one can save one random number by using only d — 1 variates to generate the diagonal 
elements. 

However, it is important to point out that an equivalent formulation is obtained by using the 
operator 



(64) 



where there is non-trivial cross correlations only between the x fluxes for Vx and y fluxes for Vy. The 
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splitting of the operator C in (64) corresponds to rewriting the the stochastic diffusion equation 
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( 62 ) in the equivalent but suggestive form 

'&t = T] V^i? + ^ V (V • 1?) 



(V • Wt) + 



Gv^v 



(65) 



where we have now distinguished between the tensorial divergence "Dt a-nd gradient operators 
Qt = T^T^ which map from tensor to vector fields and vector to tensor fields, respectively, and 
the vectorial divergence T>y and gradient operators Qy = T^v-, which map from vector to scalar 
fields and scalar to vector fields, respectively. Corresponding to the splitting of the modified 
Laplacian Cm = T^CQ = Ct + Cy into the tensorial Laplacian operator Ct = T>tQt and 



the vectorial component Cy = in Eq. (65) we have split the stochastic stress into a 

tensor white-noise field Wt in which all components are uncorrelated, and a scalar white-noise 
field Wy, which we will call the stochastic divergence stress. This representation is perhaps more 
physically-intuitive than the standard formulation in which the stochastic stress has unexpected 
exact symmetry and is exactly traceless. Note that in the more general case where the diffusion 
coefficient is spatially dependent and there is nonzero bulk viscosity rjs, the dissipative term in 



(65) becomes V- [/^(Vi?)] -|- V [{"q/S + rfs) V • with an equivalent change in the stochastic term. 
Also note that for the fluctuating incompressible Navier-Stokes equation the term with the velocity 
divergence disappears and the dissipation operator is a projected traditional Laplacian |44| [50] . 
while the stochastic flux is simply a projected tensor white-noise field. 



A. Discrete Fluctuation Dissipation Balance 

Our ultimate goal is to find a scheme that satisfies the discrete fluctuation dissipation theorem, 
that is, find a discrete modified Laplacian Lm that is a consistent approximation to the continuum 
modified Laplacian Cm{k)'d = k • C{k'd ) for small k, and a way to efficiently generate random 
increments Wm that discretize Wm and whose covariance is {WmW*^) = Lm- This task is non 
trivial in general, and completing it requires some ingenuity and insight, as illustrated in the work of 
Atzberger on multigrid methods for the scalar stochastic diffusion equation |32] . We illustrate two 
different approaches next, the first corresponding to attempting to directly discretize the modified 
Laplacian Cm, and the second corresponding to discretizing the split Laplacian Ct + Cy /'i- In the 
continuum context these are, of course, equivalent, but this is not the case in the discrete context. 
Namely, in the continuum formulation, C maps from gradients to stresses, the divergence operator 
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X> maps from fluxes to fields, and the gradient Q maps from fields to gradients. In the continuum 
context, stresses, gradients and fiuxes are all tensor fields and thus in the same Hilbert space. In 
the discrete context, however, stresses, gradients and fiuxes may be discretized differently and thus 
belong to different spaces. 



1. The modified Laplacian approach 

One approach to the problem of constructing discrete operators that satisfy the discrete 
fiuctuation-dissipation balance is to find a discretization of the divergence D and gradient G 
operators that are skew-adjoint and then form the modified Laplacian Lm = DCG = —DCD*, 
and generate the stochastic increments as Wm = DC^^^W. As discussed above, for the mean- 
ing of C^^'^ to be clear, stresses and gradients must belong to the same space. Furthermore, it 
is required that the discrete operators D and G be skew adjoint so that the discrete fiuctuation 



dissipation balance condition (30) is satisfied. 

The issue of how to define skew adjoint D and G operators also arose in the historical develop- 
ment of projection algorithms for incompressible fiow. The incompressible fiow literature suggests 
two approaches that discretize both gradients and stresses by representing them with tensors at the 
same grid of points. The first approach corresponds to fully cell-centered discretization originally 
proposed by Chorin [51 , which uses centered differences to define a skew- adjoint gradient and 
divergence operators. The second approach corresponds to a finite element-based discretization 
developed by Fortin [52] and later used in the projection algorithm of Bell et al. [53] . 

In the Fortin approach both stresses and gradients are represented as dxd tensors at the corners 
of a regular grid, where d is the spatial dimension. The divergence operator D combines the values 
of the stresses at the 2d corners of a cell to produce a value at the center of the cell. The gradient 
G = —D* combines the values of the fields at the centers of the 2d cells that share a corner into a 
gradient at that corner. In this scheme, the stochastic stresses also live at the corners of the grid. 



They are generated to have the required covariance, for example, (63). Unfortunately, the discrete 
Fortin Laplacian L = DG suffers from a serious drawback: It has a nontrivial null space. For 
example, for the scalar heat equation on a uniform grid in two dimensions, the Laplacian stencil 
obtained from the Fortin discretization is 



for which the odd {i + j odd) and even {i + j even) points on the grid are completely decoupled. In 
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Fourier space the above Laplacian is —2 [1 — cos (A/ca;) cos {Aky)] and thus vanishes for the largest 
wavevectors, \Akx\ = vr, |AA;y| = vr, which correspond to checker board zero eigenmodes. 

It can easily be verified that the same type of checker board zero eigenmodes also exist for the 
modified Fortin Laplacian Lm = DCG. In three dimensions, there are 0{N) zero eigenmodes for 
a grid of size N"^. Issues arising when using these types of stencils in the deterministic context 
are discussed in Almgren et al. Our theory for the structure factor implicitly relies on the 

definiteness of the discrete generator, and in fact, in the general non-linear setting the zero modes 
lead to instabilities of the solution of the full LLNS system of equations. We therefore abandon 
the Fortin corner-centered discretization of the fluxes. 

Fully cell-centered approximations to D and G based on second-order centered differences, 
previously studied in the context of projection methods for incompressible flows by Chorin jSlj . 
lead to a discrete Laplacian that also has a non-trivial null space and suffers similar shortcomings 
as the Fortin Laplacian. Specifically, even in one dimension one obtains a Laplacian stencil 

(L(^)n) . = [ui-2 - 2ui + Ui+2] 

where the odd-even decoupling is evident. Here we develop a cell-centered (collocated) discretiza- 
tion that preserves the null space of the continuum Laplacian. 



2. The split Laplacian approach 
An alternative to trying to form a discrete modified Laplacian Lm = Lt + Ly directly is to use 



the splitting in Eq. (65) and form the discrete tensorial Lt = DtGt and vectorial Ly = GyDy 
components separately from discretizations of the tensorial and vectorial divergence and gradient 
operators that are skew-adjoint, Gt = and Gy = Dy. The stochastic increments would simply 
be generated as DtWt + GyWy / \/?>, where Wy and the components of Wt are independent 
normal variates. 

A popular approach to discretizing the tensorial divergence and gradient operators, commonly 
referred to as a MAC discretization in projection algorithms for incompressible flow [55] , defines a 
divergence at cells centers from normal fluxes on edges, with a corresponding gradient that gives 
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normal derivatives at cell edges from cell-centered values 

(OZ)., =A.- (z<->. ^ - Z<:'. ,) + A,-' - Z<»)_ . H V ^ Z (66) 

dv 



{D*v).^i^. =Ax 1 {vi+i^j - v,j) -r — 



In this discretization, the tensor field Z 



y{x) y{x)_ y(y) ^(j/) 



is strictly divided 



into an x vector Z'^^\ which is represented on the x faces of the grid, and a y vector Z^'^\ represented 
on the y faces of the grid. The MAC discretization, which we used in the earlier one-dimensional 
examples, leads to a standard 5 point discrete Laplacian in 2D (3 point in ID, 7 point in 3D), 



whose eigenvalue in Fourier space is 2 cos {Akx) + 2 cos {Aky) — 4 and is strictly negative for all 
nonzero wavevectors, and thus does not suffer from the instabilities of the Chorin and Fortin 



discrete Laplacians, discussed in Section VII A 1 



The vectorial divergence and gradient operators cannot be discretized using the MAC frame- 
work. Namely, Dy must operate on a cell-centered vector field v, whereas the MAC-type dis- 
cretization operates on face-centered values. Instead, for the vectorial component we can use either 
the Chorin discretization [SI], in which both scalar and vector fields are both cell-centered, or 
the Fortin discretization [52j, in which scalar fields are represented at corners and vector fields 
are cell-centered. Here we choose the Fortin discretization and calculate a (scalar-valued) velocity 
divergence and the corresponding divergence stress at the corners of the grid, and also generate a 
(scalar) random divergence stress at each corner. The deterministic and random components are 
added to form the total corner-centered divergence stress, and the velocity increment is calculated 
from the (vector- valued) cell-centered gradient of the divergence stresses. Note that the nontrivial 
nuUspace of Ly does not pose a problem since Lt and thus also Lm = Lt + Ly has a trivial 
nuUspace. 

The discrete modified Laplacian that is obtained by this mixed MAC/Fortin discretization 
can be represented in terms of second-order centered-difference stencils. The first (i.e., the Vx) 
component of this Laplacian can be represented as a linear combination of the velocities in the 9 
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neighboring cells, 



r ( 



l,m=—l 
1 



^{MAC,x) (x) 



1 AMACy) (x) 

~r ~: FTJ^o — nil U 



Ay 



2 -^2-m,2+/ ^j+l,k+m, 



(67) 



r C-'^.^) (^) 



1 



where L^^^'~"'^^^^ and L^^'''^/y^ correspond to a second-order MAC and Fortin discretizations of the 
terms dxx^x and dyyi9y respectively, and L^^'^y^ discretizes dxy'&y The same stencils apply to the 
second (i.e., the Vy) component of the Laplacian as well, by symmetry, 

(J- . ^ -±- Tk''^':'y].^\ . (68) 



\y^^2-^2-m,2+'l'^j+l,k+m. + ^2-m,2+'l'^j+l,k+m 



{MAC,y)(y) 



+ 



1 j{F,xy) (x) 



3Ay2 2-m,2+ri+m,fc+i ^ 3^^^^ -^2-m,2+« ''j+m,fc+/ 

Note that we chose the peculiar indexing of the stencils so that when printed on paper they 
correspond to the usual Cartesian representation of the xy grid. The coefficients of the MAC 



stencil (67) are 



j^{MAC,x) 



while the Fortin stencils are 
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(69) 
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Results in Three Dimensions 



Our theoretical calculations have helped in formulating a complete three-stage Runge-Kutta 
scheme for solving the full LLNS system in one, two or three spatial dimensions. We have discussed 
how to generate stochastic fluxes in each stage, including the required correlations among the 
components of the stochastic stress, and have also discussed how to relate the stochastic fluxes 
in each stage. Since theoretical calculation of the three-dimensional structure factors is out of 
reach, we present some numerical results for the RK3-2RNG method in three dimensions with the 



mixed MAC/Fortin handling of the split Laplacian as given in Eqs. (67) and (68), hereafter termed 



the RK3D-2RNG algorithm. We note in passing that it is also possible to discretize the modified 
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Laplacian (see Section VII A 1 ) using a MAC-like discretization of the viscous and stochastic stresses 
that avoids the use of the Fortin corner-based discretization of the divergence stress. This saves 
one random number per ceh per stochastic flux, however, it requires the use of a non-standard 
randomized ceh-to-face projection (sphtting) of the stochastic stresses that comphcates the analysis 
and handling of physical boundaries and makes parallelization more difficult. We therefore do not 
describe this approach here, and only note that it produces very similar structure factors to those 
reported here. 

We focus on the behavior of the scheme in global equilibrium with periodic boundary conditions. 
We have implemented the full non-linear fluxes as proposed in Refs. |35| HRj . using the interpolation 



in Eq. (57) for the hyperbolic fluxes and simple interpolation of the spatially- varying viscosity and 
thermal conductivity in the handling of the viscous and stochastic fluxes. However, in the tests 
reported here we have made the magnitude of the fluctuations small compared to the means to 
ensure that the behavior is very similar to the linearized LLNS equations. Including the full non- 
linear system guarantees conservation and ensures that there are no non-linearly unstable modes. 
More careful study of the proper handling of non-linearity in the LLNS equations themselves and 
the associated numerical solvers is deferred to future publications; here, we focus on verification 
that the nonlinear scheme produces behavior consistent with the linearized analysis. We note 
that we have implemented the new RK3D algorithm also for the LLNS equations for a mixture 
of two ideal gases, closely following the original scheme described in Ref. [l8]. We find that the 
spatial discretization satisfies the discrete fluctuation-dissipation balance even in the presence of 
concentration as an additional primitive variable and that the RK3D-2RNG method performs very 
well with reasonably-large time steps. 



1. Static Structure Factors 

Examples of static structure factor Sp. for the RK3D-2RNG scheme are shown in Fig. [3j 
showing that the diagonal components S^^\ S^^'^'' and S'^'^ are close to unity, while the off-diagonal 
components S^''"'^\ g^x<'"y) ^^^^ S^'^^ are close to zero (similar results hold for S^'"'^\ not shown), 
even for a large time step (half of the stability limit). Note that the static structure factor is difficult 
to obtain accurately for the smallest wavenumbers (slowest modes) and therefore the values near 
the centers of the fc-grid should be ignored. 

It is seen in the figures that the diagonal components of are quite close to unity for the largest 
wavevectors, which is somewhat surprising, and the largest error is actually seen for intermediate 
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Figure 3: (Left) si^\ S^!'^'' and sP (top to bottom); (Right) 5, 
bottom) for RK3D-2RNG (Random Direction), with the time step a — 0.5, (3 — 3/3t/2 — 0.1, periodic 
boundary conditions with 30'^ cells, and averaging over 10® time steps. 
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wavenumbers, consistent with the one-dimensional results shown in Fig. |2] We have tested the 
method on several cell Reynolds numbers r and found that the results are worse as r increases, 
consistent with the previous analysis, however, the higher order of temporal accuracy allows for 
increasing the timestep to be a reasonable fraction of the stability limit even for large r. 

These results represent a significant improvement over the results obtained for the original RK3 
scheme presented in Bell et al. |35| HH] • Results with the original scheme were sensitive to time 
steps, requiring small time steps to obtain satisfactory results; the new scheme produces satisfactory 
results for time steps near the stability limit. Also, through the use of the mixed MAC and Fortin 
discretization, the new scheme eliminates a weak but spurious correlation S^^"^'^^^ present in the 
original scheme for small wavenumbers even in the limit of small time steps. 



2. Dynamic Structure Factors 

Examples of dynamic structure factors S^^^^ for the RK3D-2RNG scheme are shown in Fig. |4] 
as a function of lv for two relatively large wavevectors, along with the correct continuum result 
obtained by solving the system Q through a space-time Fourier transform (we did not make any 
of the usual approximations made in analytical calculations of S^ ^j [56 , and instead used Maple's 
numerical linear algebra). It is well-known that S^jf^ and S^^^ exhibit three peaks for a given k 
|56j . one central Rayleigh peak at u = similar to the peak for the diffusion equation [c.f. Eq. 



(42)], and two symmetric Brilloin peaks at a; ~ Cg/c, where Cg is the adiabatic speed of sound. 



Cs = ct-\/1 + 2/df for an ideal gas. For the velocity components, the transverse components S^"^^ 

("11 ) 

exhibit all three peaks, while the longitudinal component S^. ^ lacks the central peak, as seen in the 
figure. Note that as the fluid becomes less compressible (i.e., the speed of sound increases), there 
is an increasing separation of time-scales between the side and central spectral peaks, showing the 
familiar numerical stiffness of the full compressible Navier-Stokes equations. 

We have verified that for small wavevectors the numerical dynamic structure factors are in 
excellent agreement with the analytical predictions, even for such large time steps. For wavevectors 
that are not small compared to the discretization limits we do not expect a perfect dynamic 
structure factor, even for very small time steps. It is important, however, that the discretization 
behave reasonably for all wavevectors (e.g., there should be no spurious maxima), and be somewhat 
accurate for intermediate wavevectors, even for large time steps. As seen in Fig. [4] the RK3D- 
2RNG algorithm seems to perform well even with a large time step. Improving the accuracy at 
larger wavevectors requires using higher-order spatial differencing [57J (see discussion in Section 
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Figure 4: Diagonal (left) and the real part of the off-diagonal (right) components of the dynamic structure 
factor Sk,uj for RK3D-2RNG (dashed lines) for the same parameters as in Fig. |3] For comparison, the 
analytical solution of the LLNS equations in Fourier space are also shown (solid lines). The imaginary 
component of the off-diagonal components is less than 0.1 and it vanishes in the theory. The top part shows 
the wavevector k — (kjnax/'^j 0; 0) find the bottom shows the wavevector k — (fcmax/2, fcmax/2, k^nax/'^)- 



VC), compact stencils (linear solvers) [58], or pseudo-spectral methods [59], each of which has 
certain advantages but also significant disadvantages over the finite-volume approach in a more 
general nonlinear non-equilibrium context. 



VIII. SUMMARY AND CONCLUDING REMARKS 



In this paper we analyze finite volume schemes for the linearized Landau-Lifshitz Navier-Stokes 
(LLNS) system (j4| and related SPDEs such as the stochastic advection-diffusion equation (34 ). Our 
approach to studying the accuracy of these explicit schemes is based on evaluating the discrete static 
and dynamic structure factors, focusing on the accuracy at small wavenumber AA; = /cAx. The 



methodology for formulating the structure factor for numerical schemes is developed in sections III 



and then specialized to stochastic conservation laws in IV Applying this analysis to the stochastic 
heat equation (41) in section |v] we find the truncation error for the Euler method to be O(Atfc^); 



43 



the error for a standard predictor-corrector scheme is O(Af^fc^) using the same random numbers 
in the predictor and corrector stages but 0{At'^k^) using independent random numbers at each 



stage. Section VI extends this analysis to the third-order Runge-Kutta scheme of Beh et al. |35| 148) 
for the one- dimensional advection-diffusion SPDE. We find the best accuracy when the stochastic 



fluxes at the three stages are generated from two sets of random numbers, as given by (60); using 
this version, called RK3-2RNG, for the LLNS equations gives good results, even when nonlinear 
effects are included (see figures [2] [sj and |4| . Finally, section VII explains why the cross-correlations 
in the stress tensor in the three-dimensional LLNS require special treatment and proposes a mixed 
MAC/Fortin discretization as a way to obtain the desired discrete fluctuation-dissipation balance. 

Here we have investigated linearized PDEs with stochastic fluxes where the noise is additive. 
As such, the stability properties of the numerical schemes are the same as for the deterministic 
case. Yet in practice one would like to implement these schemes for the nonlinear stochastic PDEs 
with state-dependent stochastic fluxes. While in the limit of small fluctuations the behavior of the 
schemes is expected to be similar to the linearized case, the proper mathematical foundation and 
even formulation of the nonlinear fluctuating equations has yet to be laid out. Furthermore, the 
stability properties of numerical schemes for the nonlinear LLNS system are not well understood and 
the whole notion of stability is different than it is for deterministic schemes. For example, even at 
equilibrium, a rare fluctuation can cause a thermodynamic instability (e.g., a negative temperature 
which implies a complex sound speed) or a mechanical instability (e.g., a negative mass density). 
Capping the noises in the stochastic flux terms will not necessarily solve the problem because 
the hydrodynamic variables are time-correlated so the numerical instability may not appear on a 
single step but rather as an accumulated effect. We are investigating these issues and will discuss 
strategies to address this type of stability issue in future publications. 

One of the advantages of finite volume solvers over spectral methods is the ability to implement 
realistic, complex geometries for fluid simulations. In this paper we only consider periodic bound- 
aries but many other boundary conditions are of interest, notably, impenetrable flat hard walls 
with stick and slip conditions for the velocities and either adiabatic (zero temperature gradient) or 
thermal (constant temperature) conditions for the temperature. Equilibrium statistical mechanics 
requires that the static structure factor be oblivious to the presence of walls, even though the 
dynamic structure factors typically exhibit additional peaks due to the reflections of fluctuations 
from the boundaries. Therefore, the numerical discretization of the Laplacian operator L, the 
divergence operator D and the covariance of the stochastic fluxes C should continue to satisfy 
the discrete fluctuation-dissipation balance condition L + L* = —2DCD* and be consistent, even 
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in the presence of boundaries. Standard treatments of boundary conditions used in deterministic 
schemes can easily be implemented in the stochastic setting |35l 02] , however, satisfying the dis- 
crete fluctuation-dissipation balance is not trivial and requires modifying the stochastic fluxes and 
possibly also the finite-difference stencils near the boundaries [32]. In particular, the case of Dirich- 
let boundary conditions is more complicated, especially in the case of the mixed MAC and Fortin 
discretization of the compressible Navier-Stokes equations. Complex boundaries present further 
challenges even in the deterministic setting. We will explore the issues associated will fluctuations 
at physical boundaries in future publications. 

One motivation for the development of numerical methods for the LLNS equations is for their 
use in multi- algorithm hybrids. One emerging paradigm in the modeling and simulation of mul- 
tiscale problems is Multi- Algorithm Refinement (MAR). MAR is a general simulation approach 
that combines two or more algorithms, each of which is appropriate for a different scale regime. 
MAR schemes typically couple structurally different computational schemes such as particle-based 
molecular simulations with continuum partial differential equation (PDE) solvers. The general 
idea is to perform detailed calculations using an accurate but expensive algorithm in a small region 
(or for a short time), and couple this computation to a simpler, less expensive method applied to 
the rest. The major difficulty is in constructing hybrid is that particle and continuum methods 
treat noise in completely different ways. The challenge is to ensure that the numerical coupling 
of the particle and continuum computations is self-consistent, stable, and most importantly, does 
not adversely impact the underlying physics. These problems become particularly acute when one 
wants to accurately capture the physical fluctuations at micro and mesoscopic scales. The correct 
treatment of boundary conditions in stochastic PDE schemes is particularly difficult yet crucial in 
hybrid schemes since the coupling of the two algorithms is essentially a dynamic, two-way boundary 
condition. Recent work by Tysanner et al. [60,, Foo et al. [23 , Williams et al. 161] and Donev 
et al. have demonstrated the need to model fluctuations at the continuum level in hybrid 
continuum / particle approaches, however, a seamless coupling has yet to be developed. 

In this paper we consider the fully compressible LLNS system, for many of the phenomena of 
interest the fluid flow aspects occur at very low Mach numbers. Another topic of future work 
for stochastic PDE schemes is to construct a low Mach number fluctuating hydrodynamics algo- 
rithm. A number of researchers have considered extended versions of the incompressible Navier 
Stokes equations that include a stochastic stress tensor (TJ [251 IS] • This type of model does in- 
troduce fluctuations into the Navier Stokes equations and is applicable in some settings, such as 
in modeling simple Brownian motion. However, as pointed out by Zaitsev and Shliomis [63j, the 
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incompressible approximation introduces fictitious correlations between the velocity components 
of the fluid. Furthermore, this type of approach does not capture the full range of fluctuations in 
the compressible equations. In particular, adding a stochastic stress into the incompressible Navier 
Stokes equations creates fluctuations in velocity but does not reproduce the large scale and slow 
fluctuations in density and temperature, which persist even in the incompressible limit. We plan 
to investigate alternative formulations that can capture more of the features of the fluctuating hy- 
drodynamics while still exploiting the separation of scales inherent in low Mach number flows. We 
also note that although the theoretical importance of distinguishing between the incompressible 
approximation and the low-Mach number limit is well-established for fluctuating hydrodynamics 
\10\ |6l], numerical algorithms for the latter have yet to be developed. 
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1. Semi-Implicit Crank-Nicolson Method 

When sound is included in the fluctuating hydrodynamic equations implicit methods are not 
really beneficial since the large sound speed limits the time step. However, for the pure stochastic 
diffusion /heat equation or advection-diffusion equations with a small advection speed the time step 
may become strongly limited by the diffusive CFL limit, especially for small cells. In such cases an 
implicit method can be used to lift the diffusive stability restriction on the time step. For example, 
the second-order (in both space and time) Crank-Nicolson semi-implicit scheme for the stochastic 
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heat equation entails solving the linear system 



which is tridiagonal except at periodic boundaries. 

The analysis carried out above for explicit schemes can easily be extended to implicit methods 
since in Fourier space different wavevectors again decouple and the above iteration becomes a 
scalar linear equation for u^~^^ that can trivially be solved. Firstly, it is observed that the small 
time step limit is the same regardless of the semi-implicit treatment, specifically, the same discrete 



fluctuation-dissipation condition (30) applies. Remarkably, for the Crank-Nicolson iteration (71) it 
is found that the discrete static structure factor is independent of the time step, Sk = 1 for all /3. 



The dynamic structure factor, however, has the same spatial discretization errors (47) as for the 
Euler scheme even in the limit /? ^ 0. Furthermore, as expected, the dynamics is not accurate for 
large P and the time step cannot be enlarged much beyond the diffusive stability limit related to 
the smallest length-scale at which one wishes to correctly resolve the dynamics of the fluctuations. 

If advection is included as well also discretized semi-implicitly, the method again gives perfect 
structure factors, Sk = 1 identically, and is unconditionally stable. If only diffusion is handled 
semi-implicitly but advection is handled with a predictor-corrector approach, then it turns out 
that the optimal method is to not include a stochastic flux in the predictor step, giving the same 



leading-order error term as PC-2RNG in Eq. (53) when \r\ > 0, but giving a perfect Sk = 1 when 
r = 0. 
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